Topographic example

When we plot global data sets in 3D it can be very weird, since the data is stored with longitude-latitude angular coordinates, but the data are in metres above mean sea level. To make this work we need an arbitrary scaling to reduce the distortion due to these different units.

decimate <- function(x, dec = 10) {
  library(raster)
  r0 <- raster(x)
  res(r0) <- res(x) * dec
  setValues(r0, extract(x, coordinates(r0), method = "bilinear"))
}
library(quadmesh)
library(rgl)

## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display

## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE

scl <- function(x) (x - min(x))/diff(range(x))
data(etopo)
etopo <- decimate(etopo, 5)

## Loading required package: sp

qmtopo <- quadmesh(etopo, etopo)
open3d()

## null 
##    1

shade3d(qmtopo, col = grey(scl(qmtopo$vb[3,qmtopo$ib])), asp = c(1, 1, 0.0001))

aspect3d(1, 1, 0.1)


rglwidget()

It’s much more natural to work in a map projection, or on the surface of the sphere, but legacy spatial tools make this more difficult than it should be.

We can do what ever we like with quadmesh coordinates, without breaking the topology of the mesh at all.

library(proj4)
qmtopo$vb[1:2, ] <- t(proj4::project(t(qmtopo$vb[1:2, ]), "+proj=laea +ellps=WGS84 +lat_0=-90"))
open3d()

## null 
##    3

shade3d(qmtopo, col = grey(scl(qmtopo$vb[3,qmtopo$ib])))
aspect3d(1, 1, .1)

rglwidget()

lltopo <- quadmesh(etopo, etopo)
lltopo$vb[1:3, ] <- t(llh2xyz(t(lltopo$vb[1:3, ]), exag = 50))
shade3d(lltopo, col = grey(scl(qmtopo$vb[3,lltopo$ib])))

Try a local area.

data(etopo)
lt <- crop(etopo, extent(100, 168, -58, -40))
ltt <- ltt0 <- quadmesh(lt, lt)
ltt$vb[1:3, ] <- t(llh2xyz(t(ltt$vb[1:3, ])))
shade3d(ltt, col = grey(scl(ltt0$vb[3,ltt$ib])))

The quadmesh package also includes a mesh_plot() function for straightforward plotting of a raster to an arbitrary map projection.

plot(etopo, col = palr::bathyDeepPal(20))

prj <- "+proj=stere +lat_ts=-71 +lat_0=-90 +lon_0=147 +datum=WGS84"
mesh_plot(etopo, crs = prj, asp = 1, colfun = palr::bathyDeepPal)
xy <- proj4::project(xymap, prj)
lines(xy)